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ABSTRACT 

We introduce the N-hody simulation technique to follow structure formation in linear and nonlinear 
regimes for the extended quintessence models (scalar-tensor theories in which the scalar field has a 
self- interaction potential and behaves as dark energy), and apply it to a class of models specified 
by an inverse power-law potential and a non-minimal coupling. Our full solution of the scalar field 
perturbation confirms that, when the potential is not too nonlinear, the effects of the scalar field 
could be accurately approximated as a modification of background expansion rate plus a rescaling of 
the effective gravitational constant relevant for structure growth. For the models we consider, these 
have opposite effects, leading to a weak net effect in the linear perturbation regime. However, on 
the nonlinear scales the modified expansion rate dominates and could produce interesting signatures 
in the matter power spectrum and mass function, which might be used to improve the constraints 
on the models from cosmological data. We show that the density profiles of the dark matter halos 
are well described by the Navarro-Frenk- White formula, although the scalar field could change the 
concentration. We also derive an analytic formula for the scalar field perturbation inside halos assuming 
NFW density profile and sphericity, which agrees well with numerical results if the parameter is 
appropriately tuned. The results suggest that for the models considered, the spatial variation of 
the scalar field (and thus the locally measured gravitational constant) is very weak, and so local 
experiments could see the background variation of gravitational constant. 



1. INTRODUCTION 

The nature of the dark energy 
(Copeland, Sami & Tsujikawa 2006) is one of the 
most difficult challenges facing physicists and cosmol- 
ogists now. Although a cosmological constant (plus 
cold dark matter, to provide the concordance ACDM 
paradigm) could be a solution - and is indeed consistent 
with virtually all current observations, it suffers from 
theoretical difficulties such as why its value must be 
so small yet nonzero, and why it becomes dominant 
only at the low redshift. In all the alternative proposals 
to tackle this problem, a quintessence scalar field 
(Zlatcv, Wang & Steinhardt 1999; Wang et al. 2000) is 
perhaps the most popular one (although a new proposal 
by Barrow & Shaw (2010) provides a completely new 
type of explanation that does not require new scalar 
fields). In such models the scalar field ip is slowly rolling 
down its potential, its energy density is dominated by 
the potential energy and almost remaining constant 
provided that the potential is flat enough. The flatness 
of the potential, however, means that the mass of the 
scalar field is in general very light and as a result the 
scalar field almost does not cluster so that its effects 
in cosmology are mainly on the (modified) background 
expansion rate. 

One reason for the wide interest in quintessence 
models is that scalar fields appear in abundance in 
high-energy physics theories, in which they are often 
coupled to the curvature invariants or even other 
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matter species, leading to the so-called extended 
quintessence (Perrotta, Baccigalupi & Matarrese 

2000; Baccigalupi & Perrotta 2000; 

Baccigalupi, Perrotta & Matarrese 2000) and cou- 
pled quintessence (Amendola 2000, 2004; Jesus et al. 
2008) models respectively. The former is just a special 
class of a scalar-tensor theory (Fujii & Maeda 2003; 
Riazuelo & Uzan 2002), with the scalar field being 
the dark energy. These two classes of generalised 
quintessence models have been studied in detail in 
the linear regime in the literature (Bean 2001a,b; 
Mangano et al. 2003; Clifton, Mota & Barrow 2005; 
Nunes & Mota 2006; Pettorino, Baccigalupi & Mangano 
2005; Koivisto 2005; Brookfield et al. 2006 
Koivisto k Mota 2007; Mota & Shaw 2007, 2006 
Lee, Liu & Ng 2006; Mota et al. 2007; Bean et al. 2008 
Bean, Flanagan & Trodden 2008; Boehmer et al. 2008, 
2010). 

In recent years, studies of the cosmological be- 
haviour of the coupled quintessence model in the 
nonlinear regime have also been made, either via 
semi-analytical methods (Manera & Mota 2006; 
Mota & van dc Bruck 2004; Mota 2008; Shaw & Mota 
2008; Mota et al. 2008a; Mota, Shaw & Silk 2008; 
Saracco et al. 2009; Wintcrgerst & Pettorino 

2010), or using TV-body simulation techniques 
(Maccio et al. 2004; Nusser, Gubser & Peebles 
2005; Kesden & Kamionkowski 2006a,b; 

Springel & Farrar 2007; Farrar & Rosen 2007; 
Baldi & Pettorino 2010; Hellwing & Juszkiewicz 
2009; Keselman, Nusser & Peebles 2009, 2010; 
Hellwing, Knollmann & Knebe 2010; Baldi 2010; 
Baldi et al. 2010). In these studies the effect of the 
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scalar field is generally approximated by a Yukawa-type 
'fifth force' or by a rcscaling of the gravitational constant 
or the particle mass, without solving the scalar field 
equation explicitly. Very recently, Li & Zhao (2009, 
2010); Zhao et al. (2010); Li & Barrow (2010) gave a 
new treatment and obtained an explicit solution to the 
scalar field perturbation on a spatial grid. The new re- 
sults confirmed that the approximations adopted in the 
old literature were good for the models considered there 
(where the scalar potential was not very nonlinear), but 
for highly nonlinear potentials they broke down. 

For the extended quintessence (more generally scalar- 
tensor) models, investigations using TV-body simulations 
are rarer. The work of Pcttorino & Baccigalupi (2008), 
for example, outlined a recipe which uses certain approxi- 
mation, such as a rescaling of gravitational constant, and 
does not solve the scalar field equation of motion explic- 
itly. In Rodriguez-Meza et al. (2007); Rodriguez-Meza 
(2008a,b), the authors approximated the effect of scalar 
field coupling as a Yukawa force. However, none of these 
previous works tries to solve the scalar field on a mesh 
directly, and this is what we want to do in this work. 

The aims of this work are threefold. Firstly, we want 
to develop the formulae and methods that are needed to 
solve the scalar field explicitly, which could serve as the 
basis for future work, and to find the regime of validity 
of our method. Secondly, we want to understand whether 
the approximations adopted in the previous studies are 
good or not; given the severe limits in the computing 
power; if those approximations do work well, then one 
does not need to resort to an exact scalar field solver, 
which is considerably more economical. Finally, we want 
to study structure formation in the nonlinear regime for 
some specific models, and investigate both the scalar field 
effects on the clustering of matter and the spatial vari- 
ation of the gravitational constant (which is common to 
scalar-tensor theories). 

The organisation of this paper is as follows: In Sect. 2 
we list the basic equations which are needed in TV-body 
simulations and give their respective non-rclativistic lim- 
its. To prevent the main text from expanding too much, 
some useful expressions are listed in Appendix A, and 
the discrete versions of the resulted equations are dis- 
cussed and summarised in Appendix B. In Sect. 3 we 
briefly describe the numerical code we are using (rele- 
gating further details to Knebe, Green & Binney (2001); 
Li & Zhao (2010)), and the physical parameters of our 
simulations. We also present some results regarding the 
background cosmology and linear perturbation evolu- 
tion in our models, which could be helpful in the un- 
derstanding of the A^-body simulation results (our al- 
gorithm for the background cosmology is summarized in 
Appendix C). Sect. 4 contains the A^-body simulation re- 
sults, including key structure formation observables such 
as nonlinear matter power spectrum, mass function and 
dark matter halo profile, as well as the spatial variation 
of the scalar field. It also includes several checks of the 
approximations made in the literature. We finally sum- 
marise and conclude in Sect. 5. 

We use the unit c = 1 unless explicitly restoring c in 
the equations. The metric convention is (+,—,—,—,). 
Indices a, b, c, • • ■ run 0, 1, 2, 3 while i, j, k, ■ ■ ■ run 1,2,3. 



2. THE EQUATIONS 

This section presents the equations that will be used 
in the A-body simulations, the model parametcrisation 
and discretisation procedure for the equations. 

2.1. The Basic Equations 

We consider a general Lagrangian density for scalar- 
tensor theories 

C = -L [1 + f(<p)] R - iv>V a ^ + V(<p) - Cf, (1) 

in which = 87rG* where G* is the (bare) gravitational 
constant, R is the Ricci scalar, f(ip) is the coupling func- 
tion between the scalar field tp and curvature, V(<p) the 
potential for (p and Cf the Lagrangian density for fluid 
matter (baryons, photons, neutrinos and cold dark mat- 
ter). Note that G* is a fundamental constant of the the- 
ory. 

Varying the associated action with respect to metric 
9ab yields the energy-momentum tensor of the theory 
(note the tilde, which is used to distinguish it from the 
T a b defined below): 



T„h = T. 



f 



V a VtV5 - -^.ga&VVVc^ + g a bV((p) 



1 



— [f(ip)G ab + (ffabV C V c - VaVfc) f(tp)] (2) 



where G a b = R a b — ^g a bR is the Einstein tensor, and 
T^ b is the energy-momentum tensor for matter (including 
baryons, dark matter, neutrinos and photons, which we 
collectively refer to as 'fluid matter', although in A^-body 
simulations we use discrete particles rather than a fluid) . 
As usual, we can rearrange the Einstein equation as 



G a b '■ 

so that it now looks like 



■ Ki-kTn 



(3) 



Gab — 



1 



T f _ 
± ab 



1 



1 + / 



1 



Va<pVb<fi 



(5a&V c V c - V a V 6 ) / 
^.9ab (V<^) 2 + g a bV 



(4) 



Note the difference between T^ b and T a b\ throughout this 

paper, we will use a superscript * for normal fluid mat- 
ter, and quantities without a superscript ' always mean 
the total effective ones [the final line of Eq. (4)]. It is 
sometimes useful to define an effective Newton constant 
K e ff = re*/(l + /). Neither nor K e ff is the gravi- 
tational constant measured in a Cavendish-type experi- 
ment, which we denote instead by K0 and is given by 



2 + 2/ + 4 



df 



1 + f 



2 + 2/ + 3 



df 



(5) 



where added to make ^/kZ^p dimensionlcss, which 

is the convention we shall always follow below, it- 
self is obviously not a constant and we measure only it 
present-day value, K0 • 
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Varying the action with respect to the scalar field, tp, 
gives the scalar field equation of motion 



dVjcp) | R df(<p) 
dp 2k* dip 



:(). 



(6) 



Since wc will follow the motions of dark matter parti- 
cles in the iV-body simulations, so we also need their 
geodesic equation. The dark-matter Lagrangian for a 
point particle with mass toq is 



C 



m 



CDM 



(y) = — ^^(y - x ) A/Sabigio 



(7) 



where y is the general coordinate and Xo is the coordinate 
of the centre of the particle. From this equation we derive 
the corresponding energy-momentum tensor: 



Tiab 



m 



CDM 



<5(y-x )±g±o. 



(8) 



Taking the conservation equation for dark matter parti- 
cles (which, unlike in (Li & Zhao 2009, 2010), does not 
couple to any other matter species, including the scalar 
field tp), the geodesic equation follows as usual: 



Xq + T^ c x Xq —0, 



(9) 



where the second term on the left-hand side accounts for 
gravity. 

Eqs. (4, 6 , 9) contain all the physics needed for the fol- 
lowing analysis, though certain approximations and sim- 
plifications might have to be made in due course to make 
direct connection to A^-body simulations. 

We will consider an inverse power-law potential for the 
scalar field, 



V(tp) = 



A 4 



K*tp) 



(10) 



where a is a dimensionless constant and A is a con- 
stant with dimensions of mass. This potential has also 
been adopted in various background or linear pertur- 
bation studies of scalar fields (either minimally or non- 
minimally coupled); the tracking behaviour its produces 
makes it a good dark energy candidate and for that pur- 
pose we shall choose a ~ O(0.1 — 1). Meanwhile, the 
coupling between the scalar field and the curvature ten- 
sor is chosen to be a non-minimal one: 



f(ip)= ~fK*p 2 



(11) 



where 7 is another dimensionless constant characterising 
the strength of the coupling. Note that here again is 
added into f(tp) and V{tp) to make a dimensionless quan- 
tity ^fnZp. Although the exact value of is unknown, so 
is ip and we can solve for yj~k~lp instead of ip, not caring 
about the exact individual values of ^J~k~l and tp. 



2.2. The Non-Relativistic Limits 

The A^-body simulation only probes the motion of par- 
ticles at late times, and we are not interested in extreme 
conditions such as black hole formation and evolution, so 
we can take the non-relativistic limit of the above equa- 
tions as a good approximation. 

The existence of the scalar field and its coupling to the 
curvature leads to several possible changes with respect 
to the ACDM paradigm: 



1. The scalar field has its own energy- momentum ten- 
sor, which could change the source term of the Pois- 
son equation because the scalar field, unlike the cos- 
mological constant, can cluster (though the cluster- 
ing is often quite weak in scalar field models). Also, 
unlike in coupled scalar field models, here the X7 2 tp 
term will appear in the Poisson equation. 

2. The background cosmic expansion rate is in general 
modified, and can either slow down or speed up the 
rate of structure formation. 

3. The two gravitational potentials in the conformal 
Newtonian gauge metric ds 2 = a 2 (l + 2<p)dT 2 + 
a 2 (l — 2ip)Sijdx l dx : > , in which r and x % are respec- 
tively the conformal time and comoving coordinate, 
are no longer equal to each other (as in general rela- 
tivity), but are instead related by V 2 p (see below). 

It therefore becomes clear that the following two equa- 
tions, in their non-relativistic forms, need to be solved in 
order to obtain the gravitational force on particles: 

1. The scalar field equation of motion, which is used 
to compute explicitly the value of the scalar field tp 
at any given time and position; 

2. The Poisson equation, which is used to determine 
the gravitational potential and force at any given 
time and position from the local energy density and 
pressure, which includes the contribution from the 
scalar field (obtained from the tp equation of mo- 
tion). 

Note that unlike in the coupled scalar field models, there 
is no fifth force because there is no direct coupling to the 
particles. The scalar coupling to the curvature, however, 
does modify the gravitational potential so that gravity 
no longer follows Einstein's prescription and so this is a 
modified gravity theory. 

Wc now describe these two equations in turn. For the 
scalar field equation of motion, we denote by tp the back- 
ground value of tp and write Sip = tp — tp. Then using the 
expressions given in Appendix A we write 



a 2 \7 a \7 a tp = tp' : 



2Utp' - 
' + 3V/ 



- 2<jxp" 

- m<\>) p 1 



(12) 



in which ' = d/dr with r the conformal time, V x is 
the derivative with respect to the comoving coordinate 
x, and W = a' /a. Then, with the background part sub- 
tracted, Eq. (6) can be rewritten as 

6tp" + 2-HStp' + V 2 Jp + [V, v (tp) - V, v {0)\ a 2 

-24>(p" - (0' + 3V>' + m<j>) v>' 

[Rf v (tp)-Rf v (tp)] a 2 - 



1 



= 0, 



in which a bar denotes the background value, and the 
subscript v denotes derivatives with respect to p. Note 

that has the same sign as V^. 

In our A^-body simulations we shall work in the quasi- 
static limit, i.e., we assume that the spatial gradients are 
much greater than the time derivatives, |V x <y3| 3> |^r|- 
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Therefore, the time derivatives in the above equation are 
dropped and we obtain the simplified version 



c 2 d 2 (a<5y>) 



(13) 



\y v (<p) - v v m + ^- [rm<p) - Rum] a 3 , 



in which d 2 = — V 2 = 
convention (4-, — , — , — 



(d 2 4- dy + d% ) due to our sign 
and we have restored the factor 
c 2 in front of V 2 (the tp here and in the remaining of this 
paper is c~ 2 times the ip in the original Lagrangian unless 
otherwise stated). Note that here V has the dimension 
of mass density rather than energy density. 

To complete Eq. (13), we still need expressions for R 
and R, which are again obtained using the quantities in 
Appendix A: 



R = 

R = 
and so 



a' a 

l r 

'a 5 

a 2 a 



Qip" 4- m (0' + 3tp') - 4<9 2 ?/> 4- 2<9 2 



(14) 
(15) 



Rftp — Rftp = fipbR + RSfa 



l 



- 2 -5U (16) 
er a 



where we have again dropped time derivatives of <f> and 
ip since they are small compared with the corresponding 
spatial gradients, and SR= R— i?, Sf v = f v — f v . 

Since only <f> but not ip appears in the Poisson equa- 
tion (shown below) , we also want to eliminate the ip in 
the scalar field equation of motion. This is easy in gen- 
eral relativity, because there we have the simple relation 
4> = ip, which unfortunately no longer holds in scalar- 
tensor theories. However, we could use the i — j com- 



ponents of the Einstein equation G l 



to get a new relation between and ip. Noting that our 
A-body simulations probe the very late time evolution 
(when radiation is negligible) when the only significant 
source for (i ^ j) is the scalar field, and 



1 



(17) 



where di = d/ dx l , we could write the i — j component of 
Einstein equation as 

c 2 

did, ((f> - ip) = -Y^-jdidjf 
which gives approximately 



W-V0=-— dif 



and so 



4^-2d 2 0-2^0 4- — -ctf 

= 2^ + i^_c 2 ^. (18) 



It is important to note that in the second line of Eq. (18) 
we have implicitly linearised the equation; this is valid 
only if f(<p) is not strongly nonlinear and \8ip/ip\ <C 1- It 
turns out that the model considered in this work satis- 
fies these criteria (/(</?) <x ip 2 ). If either V(tp) or f((f) is 
highly nonlinear, then we might have \6ip\ ~ (p; in that 
case we should not approximate f to f even in the coeffi- 
cients of the perturbation variables such as d 2 Sip here, or 

write d 2 f~ f^d^Sip. The reason for the latter stricture 
is as follows: if f(<p) is highly nonlinear, then / might 
change a lot even if tp fluctuates a little, implying that 
for the linearisation to apply on our spatial grid we need 
very small grid sizes which are impossible; moreover, it 
becomes complicated to decide which solution we should 
linearise around, as the values of / in that area which 
we look at might be very different from the background 
value /. The strategy for this situation is simple: instead 
of writing d 2 f = f v d 2 Stp, we difference f(ip) directly, 
because we know the value of f(tp) in every grid cell. 
This will ensure no linearisation error. In what follows, 
however, we shall use Eq. (18), which causes negligible 
linearisation error but simplifies the equations a lot. We 
shall also write / = / in the coefficients of perturbation 
quantities such as d^.S(p and <9 2 <I>. 

Substituting Eqs. (16, 18) into Eq. (13) and rearrang- 
ing, we complete the derivation of the scalar field equa- 
tion of motion in the weak field limit, ending up with 



«.(! + /) 



c 2 <9 2 (aStp) 



= a 3 [V v {tp) - V v {(p)] - ^d 2 $ - --aSf v (19) 
for our general Lagrangian Eq. (1) and 



c <9 X (ay/K^Sip) 



1 4- 7«»( ( 9 2 
7^(^c7 2 $ - 67 {%' + H 2 ) (ay/J^Sip) 

-aK»A 4 a 3 ' 



K*ip 



, 1+a 



K*(p) 



1+a 



(20) 



for the model specified by Eqs. (10, 11), where $ = acj>. 

Next consider the Poisson equation, which is obtained 
from the Einstein equation in the weak-field and slow- 
motion limits. Here we use the — component of the 
Ricci curvature tensor, which is given as 



i?° n 



-n 2 (1-2$ 



4-3^" 4- 3H (0' 4- </>') 4- a 2 ^ 



(21) 



using the expressions in Appendix A. According to the 
Einstein equations, 



= TT (Ptot 4- 3p T OT> 2 



(22) 



where ptot and ptot are the total energy density and 
pressure, respectively. Using these two equations and 
subtracting the background part (which is just the Ray- 
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chaudhuri equation), it is straightforward to find that 



[(ptot + 3ptot) - (ptot + 3ptot)] • (23) 



in which we have dropped terms involving time deriva- 
tives of ip, </> and "H 2 </>, because they are much smaller 
than d 2 <p in the quasi-static limit. Using the energy- 
momentum tensor expressed in Eq. (4), the above equa- 
tion can be rewritten as 



(24) 



1 



1+/ 1+/ 
V{fp) V(0) 



2(1 + /) 



c z di (a5<p) 



1 + / 1 + / 



1 



1 



Ksttp ' 



— f" 

2 J 



,1 + / l + f 
for the general Lagrangian Eq. (1) and 

<9 2 $ 



- (l + 7K»^) n m Ho 



(25) 



1 



1 



1 



c 2 d 2 (ay/K^Sip) 



K*A 4 a 3 



K*A 4 a 3 



(1 + 7K*¥> 2 ) (y f K^<p) a (1 + 7K*<y5 2 ) (^kTv)" 

[(1 + 2>^)n*(p' 2 + 37K* w"] a 
1 1 



1 + 7«*<^ 2 1 + 7ft 



-r 



for the model specified by Eqs. (10, 11). In these equa- 
tions p m is the background density for matter, 5 = 
Pm/pm, and we have used the definition of Cl m given 
in Appendix C. We have also neglected the contribu- 
tion from Sip to the total density and pressure, because 
in the quasi-static limit we have \5ip"\ <C |<? 2 (5</3| and 

( 5 ( / 3 ' 2 ^ l^'^v'l \7~Lo"ip'\ ^ d 2 6tp\ (which is confirmed by 
the W-body simulation results 1 ). 

Finally, the equation of motion of the dark matter par- 
ticles is the same as in general relativity 



2— x = 

a 



1 



rV x $ 



(26) 



in which $ is determined by the modified Poisson equa- 
tion Eq. (25). The canonical momentum conjugate to x 
is p = a 2 x so from the equation above we have 



dx p 

dp 1 -. 
dt a 



(27) 
(28) 



1 According to Eq. (20) wc have (ay/nZSip) ~ O (j^^j, im- 
plying that a^fnZ&^p ~ 0(Q), so neglecting time derivatives of 8<p 
is just like dropping time derivatives of i/> and <j>, which we have 
already done to obtain the modified Poisson equation. 



Eqs. (20, 25, 27, 28) will be used in the code to eval- 
uate the forces on the dark-matter particles and evolve 
their positions and momenta in time. But before apply- 
ing them to the code we still need to switch to code units 
(see Sect. 2.3), further simplify them and create the dis- 
crete version (see Appendix B). 

2.3. Code Units 

In our numerical simulation we use a modified version 
of MLAPM ((Knebe, Green & Binney 2001)), and we will 
have to change or add our Eqs. (20, 25, 27, 28) to it. The 
first step is to convert the quantities to the code units of 
MLAPM. Here, we briefly summarise the main results. 

The (modified) MLAPM code uses the following internal 
units (where a subscript c stands for "code"): 

x c = x/B, 
Pc = p/ (H Q B) 
t c = tHo 
$ c = $/(i/ B) 2 

Pc = P/P, 

u = ac 2 ^Sp/ (HqB) 2 , (29) 

where B denotes the comoving size of the simulation box, 
Hq is the present Hubble constant, and p is the matter 
density. In the last line the quantity u is the scalar field 
perturbation Sip expressed in terms of code units and is 
new to the MLAPM code. 

In terms of u, as well as the (dimensionlcss) back- 
ground value of the scalar field, y/nip, some relevant 
quantities are expressed in full as 



A 4 



\/k<P - 



B 2 H}. 



{ B 2 H 2 \ 

/M=i+7(y^+^M 



Vf( l p) = -a- 



U (v) =27^ (^/k<p + 



1+a ' 



B 2 H 2 



(30) 



and the background counterparts of these quantities can 
be obtained simply by setting u — (recall that u repre- 
sents the perturbed part of the scalar field) in the above 
equations. 
We also define 



A = 



kA 4 
3H 2 - 



(31) 



which will be used frequently below. 

Making discrete versions of the above equations for N- 
body simulations is then straightforward, and we refer 
the interested readers to Appendix B to the whole treat- 
ment, with which we can now proceed to do A^-body 
simulations. 

3. SIMULATION DETAILS 
3.1. The N-Body Code 
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Fig. 1. — (Color Online) The background evolution in the extended quintessence models. Upper-left panel: the fractional energy densities 
for matter (f2 m ), radiation (Q r ) and the scalar field dark energy (Hoe), as indicated besides the curves, as functions of the scale factor a 
(ao = 1 today). Upper-right panel: the scalar field equation of state w = Pde/pde as a function of a. Lower-left panel: the ratio between 
the Hubble expansion rates of the extended quintessence model and ACDM as a function of a. Lower-right panel: the a-evolution of the 
effective gravitational constant that governs the growth of matter density perturbations (Gno is its value today). In all panels the black 
solid, green dot, blue dashed, purple dot-dashed, pink dot-dot-dot-dashed curves represent respectively the results for ACDM and extended 
quintessence models with (0,7) = (0.1,-0.2), (0.1,0.2), (0.5,-0.2), (0.5,0.2). 



Some of our main modifications to the MLAPM code for 
the coupled scalar field model are: 

1. We have added a solver for the scalar field, based on 
Eq. (B7). It uses a nonlinear Gauss-Seidel scheme 
for the relaxation iteration and the same crite- 
rion for convergence as the default Poisson solver. 
But it adopts a V-cycle instead of the self-adaptive 
scheme in arranging the Gauss-Seidcl iterations. 

2. The value of u solved in this way is then used 
to calculate the total matter density, which com- 
pletes the calculation of the source term for the 
Poisson equation. The latter is then solved using 
a fast Fourier transform on the domain grids and 
self-adaptive Gauss-Seidcl iteration on refinements. 

3. The gravitational potential $ obtained in this way 
is then used to compute the force, which is used to 
displace and kick the particles. 



There are a lot of additions and modifications to ensure 
smooth interface and the newly added data structures. 
For the output, as there are multilevel grids all of which 
host particles, the composite grid is inhomogeneous and 
so we choose to output the positions and momenta of the 
particles, plus the gravity and values of $ and u at the 
positions of these particles. We also output the potential 
and scalar field values on the 128 3 domain grid. 

3.2. Physical and Simulation Parameters 

The physical parameters we use in the simulations are 
as follows: the present-day dark-energy fractional energy 
density SIde = 0.743 and fi m = Ocdm + = 0.257, 
H = 71.9 km/s/Mpc, n s = 0.963, cr 8 = 0.769. Our 
simulation box has a size of 64/i -1 Mpc, where h = 
i/o/(100 km/s/Mpc). We simulate four models, with pa- 
rameters (a, 7) = (0.1,-0.2), (0.1,0.2), (0.5,-0.2) and 
(0.5, 0.2) respectively. In all those simulations, the mass 
resolution is 1.114 x 10 9 /i _1 Mq; the particle number is 
256 3 ; the domain grid is a 128 x 128 x 128 cubic and the 
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Fig. 2. — (Color Online) The CMB (left panels) and matter power spectra (right panels) for the extended quintessence models compared 
with those of the ACDM. The upper panels are for the models with a = 0.1 while the lower panels are for those with a = 0.5. The black 
solid, green dotted and blue dashed curves represent respectively the curves for ACDM and extended quintessence with 7 = —0.2 and 
7 = 0.2. For the matter power spectra, we plot the results for two different output redshifts, z = and 49, as indicated below the curves. 



finest refined grids have 16384 cells on each side, corre- 
sponding to a force resolution of about 12/i _1 kpc. 

We also run a ACDM simulation with the same phys- 
ical parameters and initial condition (see below). 

3.3. Background and Linear Perturbation Evolution 

Since the coupling between the scalar field and the 
curvature produces a time- varying effective gravitational 
constant, and the scalar field contributes to the total 
energy-momentum tensor, we expect that cosmology in 
the extended quintessence models is generally different 
from ACDM at the background and linear perturbation 
levels. A good understanding of this will be helpful in 
our analysis of the results from ./V-body simulations, and 
this is the subject of this subsection. 

Our algorithm and formulae for the background cos- 
mology are detailed in Appendix C, and are implemented 
in MAPLE. We output the relevant quantities in a prede- 
fined time grid, which could be used (via interpolation) 
in the linear perturbation and ./V-body computations. 

Fig. I shows the time evolutions of some background 
quantities of interests. For ease of comparison we have 



chosen Q m and f2 r to be the same in all models includ- 
ing the ACDM one (for definitions of fi m and f2 r see 
Appendix C), and as a result in the upper left panel the 
curves for different models converge at common right- 
hand ends. We see increasing a results in an earlier and 
slower growth of Ode (^de = 1 — — £l r )- This indi- 
cates a larger dark energy equation of state parameter, w, 
which is confirmed by the upper right panel. Physically, 
this is because, the larger a is, the steeper the potential 
becomes and thus the faster the scalar field rolls. Notice 
that w is also larger for positive 7, with a being the same. 
This is because in Eq. (6) the Ricci scalar R < and for 
positive 7 the term ^"/v has the same sign as V v , thus 
helping the scalar field to roll faster. Because of its large 
predicted value of w, the model (a, 7) = (0.5,0.2) is al- 
ready excluded by cosmological data, but here we shall 
keep it for purely theoretical interest (i.e., to see how 
changing a or 7 changes the nonlinear structure forma- 
tion). 

We are also interested in how the expansion rate in an 
extended quintessence model differs from that in ACDM, 
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Fig. 3. — (Color Online) The relation between the magnitudes of the scalar field perturbation ay/ K„S<p (in unit of 10 7 ) and gravitational 
potential <E> (in unit of 10 — 6 ) for the four extended quintessence models (the four columns) at three different output times a = 0.3,0.5, 1.0 
(the three rows) as indicated above the frames. The black solid line in each panel represents the analytical approximation Eq. (35) (see 
text) and the ~ 10,000 green dots the results from a thin slice of our simulation boxes. 



and the results for our models are shown in the lower-left 
panel of Fig. 1 , which plots the H/Hacdm as a function 
of a. The rather odd behaviour of the models at low 
rcdshift is because of the complicated evolution of the 
scalar field (and the fact that we have chosen Ho to be 
the same for all models, again for ease of comparison), 
while the high-redshift behaviour could be seen directly 
from Eq. (C4). In Eq. (C4) the energy density of the 
scalar field can be dropped at high z, and so we have 



U_ 

Ho 



l + fo 
1 + / 



(32) 



where we have also neglected the radiation for simplicity 
(which is valid after the matter-radiation equality) . This 
shows that in extended quintessence models the gravita- 
tional constant relevant for the background cosmology is 
rescaled by (l + /o)/(l + /). Because fo = /(</?o) where ipo 
is the present-day value of ip, and (p is monotonically in- 
creasing in time, so for our choice of /((/?) [cf. Eq. (11)] we 
have (l+/ )/(l+/) > 1 for 7 > and (l+/ )/(l+/) < 1 
for 7 < 0: thus models with 7 > have H/Hacdm > 1- 
It turns out that the gravitational constant relevant for 
the growth of matter density perturbations is also differ- 



ent from the one governing the background cosmology. If 
we denote the matter density perturbation by S rn , then 
it can be shown, using the linear perturbation equations, 
that on small scales the evolution equation for 6 m reduces 
to 



HS' — G^—-^-^t m 8 m a 2 



(33) 



in which ' = d/dr and r is the conformal time (see Ap- 
pendix C), and we have defined 



Gn = 



l + /o 



2 + 2/ + 4 



AS 



dy/K* if 



1 + / 



2 + 2/ + 3 



df 



(34) 



Note that this quantity could also be directly read off 
from the modified Poisson equation Eq. (B3). 

In the lower right panel of Fig. 1 we display the evolu- 
tion for Gn in the models considered. Again, Gn is larger 
at earlier times for positive 7 and smaller for negative 7, 
because of our specific choice of f{ip) in Eq. (11), and 
the fact that ip is always increasing in time. 

It is well known that a higher rate of background ex- 
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Fig. 4. — (Color Online) The relation between the magnitudes of the naive gravity (see text) and full gravity for the four extended 
quintessence models (the four columns) at three different output times a = 0.3, 0.5, 1.0 (the three rows) as indicated above the frames. The 
black solid line in each panel represents the analytical approximation Eq. (34) (see text) and the ~ 10, 000 green dots the results from the 
simulations. 



pansion means that structures have less time to form, and 
a larger Gn speeds up the structure formation. These two 
effects therefore cancel each other to some extent, which 
results in a weaker net effect of an extended quintessence 
field on the large scale structure formation. This is con- 
firmed by our linear perturbation computation depicted 
in Fig. 2. In the right-hand panels of this figure we have 
plotted the matter power spectra for different models at 
two different redshifts (0 and 49). It is interesting to note 
that on small scales the matter power is closer to that of 
ACDM, despite the significant differences in background 
expansion rate and Gn (cf. Fig. 1). Because of this, we 
shall choose ACDM initial condition for our iV-body sim- 
ulations for all our models, saving the effort of generating 
separate initial conditions for different models. 

The left hand panels of Fig. 2 display the CMB power 
spectra for the models we consider. Again the difference 
from ACDM is fairly small, and there is only a small 
shift of the CMB peaks even though the background 
expansion rate changes quite a bit. The latter is be- 
cause peak positions are determined by the ratio of the 
sound horizon size at decoupling and the angular dis- 
tance to the decoupling, and in our model both of these 
decrease/increase as the Universe expands faster/more 
slowly, their ratio does not change much. 

To briefly summarise, the study of background cosmol- 
ogy and linear perturbation shows that a modified back- 



ground expansion rate and a rescaled gravitational con- 
stant, the two most important factors affecting structure 
formation in extended quintessence models are opposite 
effects. It is then of interest to see how these two effects 
compete in the nonlinear regime. 

4. N-BODY SIMULATION RESULTS 

This section lists the results of extended quintessence 
A^-body simulations. We shall start with a few prelimi- 
nary results which both give some basic idea about the 
extended quintessence effects and serve as a cross check 
of our codes. Then we discuss the key observables for 
the nonlinear structure formation such as matter power 
spectrum, mass function and halo properties. We also 
comment on the halo profile of the scalar field and the 
spatial variation of gravitational constant. 

4.1. Preliminary Results 

As mentioned above, in both the linear and A^-body 
codes we compute background quantities via an interpo- 
lation of some pre-computed table. Because background 
cosmology is important in determining the structure for- 
mation, it is important to check its accuracy. For this we 
have recorded in Table 1 The age of the universe today 
for different models as computed by these two codes. The 
two codes are compatible with each other indeed. 
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Fig. 5. — (Color Online) The fractional difference between the nonlinear power spectra for the extended quintessential and ACDM models. 
The results for the four models of (a, 7) = (0.1, —0.2), (0.1, 0.2), (0.5, —0.2) and (0.5, 0.2) are respectively represented by the black, green, 
purple and pink curves. The four panels are for four output times a = 0.3, 0.5, 0.7 and 1.0 as indicated inside the corresponding frames. 



TABLE 1 

Current age of the Universe for the different models 
under consideration as computed by the linear 
perturbation and tv-body codes. unit is gyr. 



model 



linear code TV-body code 



ACDM 13.680 13.678 

{a, 7) = (0.1,-0.2) 13.639 13.638 

(a, 7) = (0.1,0.2) 13.408 13.408 

(a, 7) = (0.5,-0.2) 13.513 13.513 

(a, 7) = (0.5,0.2) 12.097 12.096 

Because one of the advantages of our TV-body code is 
that it solves the scalar field perturbation explicitly, it is 
important to check that the solution is with expectations. 
From Eqs. (Bl, B2) it could be seen clearly that, if the 
contribution to the local density and pressure from the 
scalar field is negligible compared with that from matter, 
then the modified Poisson equation and scalar field equa- 
tion of motion end up with the same source term (up to 



a (^-dependent coefficient). In this situation we expect 

27^k7<^ 



1 



(35) 



which means that u is simply proportional to $ c with 
a time-dependent coefficient. In Fig. 3 we have checked 
this relation explicitly: we select a thin slice of the sim- 
ulation box, fetch the values for u and 4> c at the posi- 
tions of the particles (about 10000 in total) therein, and 
display them as scatter plots. The solid curve is the ap- 
proximation Eq. (35) while the green dots are simulation 
results; we can see they agree very well with each other, 
showing that the above approximation is a good one. 
Note that the scalar field perturbation a^/x^Sip is gener- 
ally less than 10 -6 , compared with the background value 
yfnZCp ~ O(0.1 — 1). This confirms that it is consistent to 
neglect the perturbation in scalar field density/pressure, 
drop terms such as Sip and Sip, and replace tp by <p in co- 



efficients of perturbation quantities such as d£ 



■Sip) 



and d^Q. It also serves as a check of the numerical code. 

As a final consistency check, let us consider the total 
gravitational force on particles. In extended quintessence 
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Fig. 6. — (Color Online) The mass functions for the mod- 
els considered. The black solid, green-dotted, bluc-dashed, purple 
dot-dashed and pink dot-dot-dot-dashed curves stand for the re- 
sults for ACDM and extended quintessence models with (0,7) = 
(0.1,-0.2), (0.1,0.2), (0.5,-0.2) and (0.5,0.2) respectively. The 
horizontal axis denotes the halo mass (in unit of h^ 1 Mg) and the 

vertical axis is the halo number density (in unit of /i 3 Mpc — 3 ). Only 
the results at a = 1 are plotted. 

models, this is given by Eq. (B2), and when the pertur- 
bation in the scalar field density/pressure is negligible 
(which is the case as shown above) we get 

V 2 $ c «^G N r> m ffo (Pc-l) (36) 

in which Gn is given in Eq. (34). On the other hand, if 
we consider (naively) that gravity is described by general 
relativity, then we should neglect the Gn on the right- 
hand side. Manipulating Eqs. (Bl, B2) we obtain: 

\±^4v> ( * c + j^^u) « \n m Hl iPc - 1) . 

l + -fK*ip z \ l+7ft*<^ / 2 

Thus (* c + i+^Zt' 2 u ) acts as thc P° tcntial for 

naive gravity (i.e., general relativity), and by differcing 
it we could obtain the naive gravitational force. In Fig. 4 
we show the scatter plot of the naive gravity versus full 
gravity for the same particles as in Fig. 3 (green dots) 
as well as their approximate ratio Gn (solid line). Again, 
the agreement is remarkably good. 

4.2. Nonlinear Matter Power Spectrum 

As we have seen above, the linear matter power spec- 
trum for the extended quintessence model really does 
not show much useful information on small scales, and 
so we need to investigate whether nonlinear effects could 
change this situation and therefore potentially place 
more meaningful constraints. 

Fig. 5 provides a positive answer to this question. 
Here we have plotted the fractional difference of the ex- 
tended quintessential nonlinear matter power spectrum 
from that for ACDM (remember that we use the same 
initial condition for all simulations). We can see that for 
the models with a = 0.1 thc differences are small even in 
the nonlinear regime, indicating that the scalar field re- 
ally does not affect the matter distribution significantly 
if the potential is flat. However, for thc a = 0.5 cases 



in which the coupling strength 7 remains the same, the 
difference could be as large as 30% ~ 50%, guaranteeing 
an observable signature. 

Furthermore, for negative 7 (the purple curve) the ex- 
tended quintessential power spectrum beats the ACDM 
one on small scales, whereas for the positive 7 case (the 
pink curve) it is just the opposite. As shown before, when 
7 < 0, both the background expansion rate and the ef- 
fective gravitational constant governing the structure for- 
mation decrease, boosting and weakening the collapse of 
matter respectively. In our a = 0.5 cases the first effect 
has clearly taken over on small scales. 

4.3. Mass Function 

A second important observable is the mass function. 
This gives the number density of dark matter halos as a 
function of halo mass. For this we need to identify thc 
dark matter halos from thc output particle distribution 
of the A-body simulations, and this determination is per- 
formed using a modified version of MHF (Knebe & Gibson 
2004), MLAPM's default halo finder. 

MHF optimally utilizes the refinement structure of thc 
simulation grids to pin down the regions in which po- 
tential halos reside and organize the refinement hierar- 
chy into a tree structure. MLAPM refines grids according to 
the particle density on them and so the boundaries of the 
refinements are simply isodensity contours. MHF collects 
the particles within these isodensity contours (as well 
as some particles outside). It then performs thc follow- 
ing operations: (i) assuming spherical symmetry of thc 
halo, calculate the escape velocity v esc at the position of 
each particle, (ii) if the velocity of the particle exceeds 
v esc then it does not belong to the virialized halo and is 
removed. Steps (i) and (ii) are then iterated until all un- 
bound particles are removed from the halo or the number 
of particles in the halo falls below a pre-defined threshold, 
which is 20 in our simulations. Note that thc removal of 
unbound particles is not used in some halo finders using 
the spherical overdensity (SO) algorithm, which includes 
the particles in the halo as long as they are within thc 
radius of a virial density contrast. Another advantage of 
MHF is that it does not require a predefined linking length 
in finding halos, such as the friend-of-friend procedure. 

Our modification to MHF is simple: because the effective 
gravitational constant in the extended quintessence mod- 
els is rescaled by a factor Gn [cf. Eq. (34)], the escape 
velocity of particles from a halo is also multiplied by this 
factor, and in MHF we have only changed the criterion for 
removing particles from virialised halos accordingly. In 
reality, because we are only interested in the a = 1 halos 
in this work, Gn is quite close to 1 and the effect of our 
modification is not large. 

The mass functions for our simulated models arc shown 
in Fig. G. It shows that all extended quintessence models 
considered here, irrespective of their parameters, produce 
less massive halos than ACDM, whereas (only) the model 
(a, 7) = (0.5, —0.2) produces a larger number of less mas- 
sive halos. These features are in broad agreement with 
those shown in the matter power spectra (Fig. 5) where 
all models show less matter clustering on the large scales, 
whereas (only) the model (a, 7) = (0.5, —0.2) shows more 
power on small scales. Thc physical reason is again the 
competition between the modified background expansion 
rate and rescaled effective gravitational constant Gn- 
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Fig. 7. — (Color Online) The NFW fitting results for two halos randomly selected from the 80 most massive halos in each simulation (see 
text for details). The upper and lower green asterisks represent respectively the density profile from TV-body simulation for the more and 
less massive halo, and the green solid curves their NFW fittings. For comparison we also shown the corresponding TV-body (black crosses) 
and fitting (black dashed curves) results for the ACDM model. The horizontal axis is the distance from halo centre (in units of Ti~ x kpc) 
and the vertical axis is the density contrast. The four panels are for the four models as indicated above the frames. 



4.4. Halo Properties 

In the ACDM paradigm, it is well known that 
the internal density profiles of dark matter halos 
are very well described by the Navarro-Frenk- White 
(Navarro, Frcnk & White 1996) formalism 



D~( \ 2 ^ ' 

where p c is the critical density for matter, /3 is a dimcn- 
sionlcss fitting parameter and R s & second fitting param- 
eter with length dimension. j3 and R s are generally differ- 
ent for different halos and should be fitted for individual 
halos, but the formula Eq. (37) is quite universal. 

We are thus interested in whether the halo profiles in 
an extended quintessential Universe are also featured by 
this universal form. For this we select the 80 most massive 
halos from each simulation and fit their density profiles 
to Eq. (37). The results show that the NFW profile de- 
scribes the extended quintessential halos at least as well 
as it does for the ACDM halos. Fig. 7 shows the fittings 
for two halos randomly picked out of the 80: one at ~ 
(10.34, 28.63, 13.91)/i -1 Mpc with mass - 1.88 x 10 14 M Q 



and the other at (41.77, 31.91, 21.20)/i _1 Mpc with a mass 
- 4.98 x 10 13 M O - 

There are some interesting features in Fig. 7. Firstly, 
for the models with a = 0.1 (the top panels) the halo 
density profile for extended quintessence models (green 
asterisks) is very similar to the ACDM results (black 
crosses) and thus their fittings almost coincide. Secondly, 
for the model of (a, 7) = (0.5,-0.2), the chosen ha- 
los show more concentration of the density profiles in 
the scalar model than in ACDM. Thirdly, the model of 
(a, 7) = (0.5, —0.2) has just the opposite trend and suf- 
fers a suppression of density in large parts of chosen ha- 
los. 

To verify that the above features are actually typical 
for the corresponding models, we have plotted in Fig. 8 
the fitting results for all the 80 massive halos in all simu- 
lated models. Here in addition to the NFW concentration 
parameter cnfw = r 200 /Rs, where r2oo is the radius at 
which the density is equal to 200 times the critical den- 
sity p c and R s the NFW parameter, we have also shown 
the fitting errors for each halo. 

We would like to point out several important impli- 
cations of Fig. 8. Firstly, for all models the fitting error 
for the extended quintessential halos (lower green aster- 
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Fig. 8. — (Color Online) Scatter plot of the NFW fitting of the dark matter halo density profiles for the 80 most massive halos in each 
simulation box. In all panels the black crosses and green asterisks stand for results of ACDM and extended quintessence models respectively. 
The upper cluster of points in each panel represents the fitted cnfw and the lower cluster is for the fitting error, as indicated beside the 
vertical axis. The horizonal axis is the halo mass (in units of /i — 1 Mg). The four panels are for the models as indicated above the frames. 



isks) is comparable to that for the ACDM halos (lower 
black crosses), indicating that the density profiles for the 
former are equally well described by the NFW formula 
Eq. (37). Secondly, for the models with a = 0.1 (the top 
panels) we can see that the fitted cnfw for the extended 
quintessential halos is comparable to that for ACDM, 
which is in agreement with our finding in Fig. 7 that the 
density profiles for the chosen halos are almost the same 
as in the ACDM prediction. Thirdly, for the model of 
(a, 7) = (0.5,-0.2), the halos tend to be more concen- 
trated (i.e., with larger cnfw) than in ACDM. Fourthly, 
for the model (a, 7) = (0.5, 0.2), the halos tend to be less 
concentrated (i.e., with smaller cnfw) than in ACDM. 
The above three features show that our qualitative find- 
ings in Fig. 7 are quite typical. Finally, the halo masses in 
the model of (a, 7) = (0.5, —0.2) are on average smaller 
than those in ACDM, because the upper green asterisks 
in the lower right panel consistently shift leftwards with 
respect to the upper black crosses: this is consistent with 
the mass function result that this model produces less 
massive halos than ACDM. 

In summary, the halo density profiles for the extended 
quintessence models are well described by the NFW for- 
mula, but the existence of the scalar field and in particu- 



lar its coupling to curvature do change the concentration 
parameters of the halos, so long as the potential is not 
too flat. It seems that the modified background expansion 
rate beats the effect of the rescaled effective gravitational 
constant here. 

4.5. Halo Profile for Scalar Field Perturbation 

We have already seen that the coupling between the 
scalar field and the curvature scalar causes time and spa- 
tial variations of the locally measured gravitational con- 
stant K0 . It is then of our interest to ask how Km varies 
across a given halo and whether this could produce ob- 
servable effects. This subsection answers this question, 
by giving an analytical formula and comparing it with 
numerical results. 

Recall that Fig. 3 shows that to a high precision the 
scalar field perturbation a^/n^Sif is proportional to the 
gravitational potential $ [cf. Eq. (35)] everywhere. This 
means that if we could derive an analytical formula for <fr 
in halos, then we know ay/~K^6(p straightforwardly. Such 
a derivation has been done in (Li, Mota & Barrow 2010) 
for a different model, but here we shall briefly repeat it 
for the extended quintessence model for completeness. 

Assuming Eq. (37) as the density profile and sphericity 
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Fig. 9. — (Color Online) Comparison between our analytic formula for the scalar field perturbation a^fnZ&ip [Eq. (43)] and the results 
from numerical simulation. The nine panels are for nine halos selected from the simulation box, whose masses are indicated inside each 
frame. In each panel the solid curve is Eq. (43) with <E>* = 0, green crosses are the numerical results for dy/n^&tp, the dashed curve is 
Eq. (43) with <3>* appropriately tuned to match the green crosses, and the red asterisks are the a-^/reTtSi/j computed from the value of <f> 
using Eq. (35). The horizontal axis is the distance from the halo centre, and vertical axis stands for the value of a^/KZ&^p. 



of halos, we can derive V c (r), the circular velocity of a 
particle moving around the halo at a distance r from halo 
centre, to be 



GM(r) 



iln 

r 



1 + i 



R s 



(38) 



where M(r) is the mass enclosed in radius r, G is 
the properly rescaled gravitational constant. Again, this 
equation is parameterized by (3 and R s . From a simula- 
tion point of view, it is straightforward to measure M(r) 
and then use Eq. (38), instead of Eq. (37), to fit the val- 
ues of (3 and R s ; from an observational viewpoint, it is 
easy to measure V c (r), which could again be used to fit 
(3 and R s . 

The potential inside a spherical halo is then given as 



$(r): 



GM(r'] 

Z?2 



dr' + C 



(39) 



in which GM(r)/r 2 is the gravitational force and C is a 
constant to be fixed using the fact that <I>(r = oo) = $00 
where 3>oo is the value of the potential far from the halo. 
Using the formula for GM(r)/r 2 given in Eq. (38) it is 



not difficult to find that 



' "^Pdr< =A,G PPc Rl 



hl ( 1 + 7C 



R, 



and so 

C=$ 00 -4nGpp c R 2 s . 
Then it follows that 

$(r) = $co - ^Gf3 Pc ^ In ( 1 + -J- 
r V R s 



If the halo is isolated, then $ 
$(r) 



and we get 



R 3 ( r 
i7rGj3 p c ^Ll n l 1 + 

r V R s 



(40) 
(41) 

(42) 



However, in iV-body simulations, we have a large number 
of dark matter halos and no halo is totally isolated from 
the others. In such situations, in Eq. (41) should be 
replaced by $*, which is the potential produced by other 
halos inside the considered halo (note that in practice $* 
could be position dependent as well, but for simplicity we 
assume that it is a constant, which is a good assumption 
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for many halos). Then we get 



1 



R 3 ( T 

-47rG/?p c ^ln 1 + — 
r V R s 



(43) 



Eq. (43) provides a neat analytical formula for a^/TiStp 
in halos, but unfortunately in most cases it cannot be 
used directly because we lack information about $*. We 
will then be forced either to fit <£>„, as a free parameter, 
or tune its value to match simulations or observations. 
In this work we shall take the second approach, and we 
find that with an appropriate value of <I>* and with values 
of (3 and R s fitted using Eq. (38), Eq. (43) agrees with 
numerical results for most halos. 

Some examples are shown in Fig. 9, in which we have 
computed a^/K,5tp(r) using four different methods: direct 
TV-body simulation results (big green crosses), Eq. (43) 
with $» = (solid curves), Eq. (43) with properly 
tuned (dashed curves) and Eq. (35) with $ directly from 
.ZV-body simulations (small red asterisks). Clearly the 
crosses and asterisks agree with each other very well, 
which is another demonstration that Eq. (35) is a very 
good approximation (cf. Fig. 3). The solid curves differ- 
ent significantly from the numerical results, showing that 
$* is actually nonzero; once it is appropriately tuned, 
then Eq. (43) (dashed curves) agree with the numerical 
results very well for all the chosen halos. Eq. (43) there- 
fore provides a useful analytical formula which might aid 
in general analysis. 

We also notice that across the halos, the variation of 
ay/^Stp is typically < O (l(T 6 - 1(H). Such a small 
variation is unlikely to be detectable using current obser- 
vational instruments, and thus we do not expect special 
constraints based on the spatial variation of G. How- 
ever, we stress that the above result is only for a class 
of extended quintessence models, and although we ex- 
pect it to be valid for other potentials which are not 
particularly nonlinear, the situation could be dramati- 
cally changed in cases where the potential or coupling 
function becomes highly nonlinear. Such models require 
a more careful treatment, including some of the approx- 
imations adopted above becoming invalid, and are thus 
beyond the scope of the current work. 

5. SUMMARY AND CONCLUSION 

In summary, in this paper we have described a nu- 
merical method to study extended quintessence models, 
where the quintessence field has a scalar-tensor type of 
coupling to the curvature, from background cosmology to 
nonlinear structure formation, and discussed the regime 
of validity of the method. Instead of assuming a Yukawa 
force due to scalar coupling or simply a rescaling of grav- 
itational constant, we have solved the scalar field and its 
spatial variation explicitly from their equation of motion. 
This is a necessary step in general to obtain trustable re- 
sults and check various approximations which are made 
to simplify the computation. 

As specific examples, we apply the above method to 
a specific class of models with inverse power-law poten- 
tial Eq. (10) and non-minimal coupling Eq. (11). The 
analysis of the background cosmology and its linear per- 



turbation shows that for these models the effective grav- 
itational 'constants' relevant for the cosmic expansion 
rate and structure formation are either both increased 
or both decreased (albeit by slightly different amounts). 
The two effects compete and cancel each other, and as a 
result the net effect on large scale structure in the linear 
regime is weak (cf. Fig. 2). We then investigated whether 
a more significant signature of the scalar field could be 
imprinted in the nonlinear regime of structure formation. 

The nonlinear matter power spectra plotted in Fig. 5 
suggests that the effect of the scalar field is more sig- 
nificant in the nonlinear regime. For the models with 
a = 0.5 (i.e., steeper potential), the scalar field changes 
(cither increases or decreases) the matter power spec- 
trum by 30 ~ 50% on small scales with respect to the 
ACDM prediction. Going to nonlinear scales thus greatly 
enhances the power of constraining such models using 
cosmological data. However, the power is more limited 
for models with a = 0.1 (i.e., shallower potential); their 
matter power spectra are very similar to the ACDM re- 
sults. Of the two competing effects mentioned above, we 
find that the modified background expansion rate is more 
influential on nonlinear scales. 

Properties of mass functions (cf. Fig. G) are in quali- 
tative agreement with what we have seen in the matter 
power spectrum, with the extended quintessence mod- 
els producing less massive halos than ACDM. Therefore 
galaxy cluster counts could place meaningful constraints 
on such models as well. But as the matter power spec- 
trum, the mass function for the models with a = 0.1 
(i.e., shallower potential) is very similar to the ACDM 
result. 

The halo density profiles for the extended quintessence 
models are shown to be well described by the well-known 
NFW formula (cf. Fig. 7, 8). In Fig. 8 we have shown the 
results of the fitting for the 80 most massive halos from 
each simulation. Consistent with the findings in Figs. 5 
and 6, we see that the concentration parameter Cnfw for 
the halos in the a = 0.1 models is almost the same as 
for the ACDM halos. But for a = 0.5, the 7 = -0.2 and 
7 = 0.2 cases predict overall bigger and smaller cnfw 
than ACDM, respectively. Furthermore, Fig. 8 shows 
clearly that the halos in the (a, 7 ) = (0.5,0.2) model 
are consistently less massive than those in ACDM, as 
suggested by the mass function plots. 

Scalar-tensor theories (which the extended 
quintessence models belong to) are often studied 
in the context of varying gravitational constant, and so 
we have also considered the spatial variations (time vari- 
ation has been investigated in detail elsewhere and will 
not be repeated here) in the scalar field (or equivalcntly 
the locally measured gravitational constant k^). We 
first showed in Fig. 3 that the approximation that the 
scalar field perturbation a^/n^Sip is proportional to the 
gravitational potential $ [cf. Eq. (35)] is fairly accurate. 
Then, based on this fact and using the NFW density 
profile, we derive an analytical formula for a-JK^8(p(r) in 
spherical halos, in which the parameters are obtained by 
fitting the NFW circular velocity profile. We have shown 
that this formula could be tuned to fit the numerical 
results pretty well for most halos (cf. Fig. 9). 

Fig. 9 indicates that the spatial variation of a^/Ti^Stp 
across halos is at most of order 10" 5 . which is far smaller 
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than the background value y^ip ~ 0(0.1 — 1). Therefore 
the spatial variation of Km is expected to be of order 10 -5 
or less in the halos, which is difficult to detect. 

The smallness of a^/n^Sif also implies that the approx- 
imations we have made to simplify the simulations are 
valid. For example, because | a^fRZS tp\ <C 1, which means 
it is reasonable to ignore the contribution from Sip, Sip to 
the total density/pressure perturbation, we can also re- 
place ip by if in the coefficients of perturbation quantities 

such as d 2 $> and d 2 [a^/n^Sip). Moreover, the quasi static 

limit, i.e., neglecting 8ip,5f compared to d 2 {a^J^Sip) , 
is guaranteed to work well. 

One of the most important results of this work is 
that it confirms explicitly that, for a broad range of ex- 
tended quintessence models, the TV-body simulation re- 
duces to modifying the background expansion rate and 
rescaling the effective gravitational constant based on the 
the background value of ip. This works to quite high 
accuracy and thus there is no need to solve the scalar 
field equation of motion explicitly, which is particularly 
time-consuming for large simulations. However, we ex- 
pect this approximation to break down in extreme situa- 



tions where the potential (or perhaps the coupling func- 
tion) becomes highly nonlinear, and then both our results 
and method might have to be revised. 
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power spectrum from output particle distribution, and 
a modified version of CAMB (Lewis, Challinor & Lascnby 
2000) for our linear perturbation computation. We thank 
David Wands for discussions. B. Li is supported by 
the Research Fellowship at Queens' College, Cambridge, 
and the Science and Technology Facility Council of the 
United Kingdom. DFM thanks the Research Council of 
Norway FRINAT grant 197251/V30. 



APPENDIX 



USEFUL EXPRESSIONS 

In this appendix we list some useful expressions in the derivation of our equations, because different researchers use 
different conventions. 
Our line element is 



ds 2 =a 2 (l + 2(f>)dT 2 - a 2 (l - 2tp)-f ij dx i dx j 



(Al) 



where r is the conformal time, x l is the comoving coordinate and 7^ the metric in the 3-space (with i, j running over 
1,2,3). The nonzero Christoflc symbols, up to first order in perturbation, are 



1 00 — 



4>\ 



0i 



0j 



r?. = -(l-20-2VO7i,--V' , 7y 



(A2) 



where a comma denotes a partial derivative with respect to the comoving coordinate, and indices are raised and lowered 
by 7 1 - 7 and 7^ respectively. ' = d/dr. 
The Ricci tensor is 



ti-ab — 1 ab.c ~ 1 ac.b "+" 1 cd 1 ab ~ 1 cb L ad 



^d tic 



(A3) 



and its components up to first order in perturbation are 



R< 



no : 



a 
a 



+ + 3- {(/)' +i>'). 
a 



^ = 2^ + 2-^ 
a 



-i>"n j --{4>' ' + 5^')7« -{4>- VO, 



a 
a 



(l-2<£-2V)7y +t k lir 
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The Ricci scalar R and relevant components ol Einstein tensor G a b = R a b — hgabR are 



R 



6 



_£_(!_ 20) + 1(^ + 3^0 
a a 



a z \ a / a z a or ' 



-2V/'-2-(0' + 2V/)-(0-# fe fe 



(1 - 20)<5* 



(A4) 
(A5) 



(A6) 



DISCRETE EQUATIONS FOR THE iV-BODY SIMULATIONS 

In the MLAPM code the Poisson equation Eq. (25) is (and in our modified code the scalar field equation of motion 
Eq. (20) will also be) solved on discretised grid points, so we must develop the discrete versions of Eqs. (20, 25) to be 
implemented in the code. Before doing that, we note that Eqs. (20, 25) are not independent but are coupled together, 

which could further complicate the solver. As a result, we should first decouple them by eliminating d 2 (a,y/K^5(p) 

(d 2 &) from the equation for $ (Sip). This is easy to do and the resulted equations are respectively 



2^ ,7:2 i 



67 K^ip 



1 + 7K*(/3 



c 2 di (ay/J^Sip) = -67 {%' + H 2 ) ay^Sip - 3aXH'^a 3 



1 



Sjy/K^ip (l + 7K*<y5o) ^m#0 



P< 



1 



+67 y/K^(fi XHq a 3 



1 + 7K*!y2 2 1+7K*<£ 2 

1 



(1 + "fK*lf 2 ) (i/K^ip) a (1 + 7K»y 2 ) (^/K^ifY 

1 1 



-2^^/nlpa [(1 + 37)k*<// 2 + S^n^pip' 



(Bl) 



for the scalar field, and 



1 + JK*p 2 + 6-f 2 K*,<fi 2 



1 + "/K*p 2 + 8j 2 K*p 2 X 



- (1 + jn*(p ) n m H 



Pc 



1 



+a [(1 + 37)k*</3 /2 + 3jK^ipip" 

' -1 — 9 o 9 — 9 0^ 

1 + 7«*v 3 + °7 K*cp 



1 + 7«*(/? 2 1 + JK^tp 2 

1 



2„3 



3XH£a 



1 



l+7K*ly9 2 1+7K*C/? 2 
1 1 



(1 + "fK t p 2 ) (y/K^p) a (I + JK^p 2 ) (i/K^ip) 1 
6j 2 y/K^p 



j (W + H 2 ) ay/^Sip 
1 + 7K*<^ 2 + 87 2 k*^ 2 v ' v 



*0 J 



(B2) 



for the gravitational potential. 



Introducing the variable u (cf. Sect. 2.3), the Poisson 
equation becomes 
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1 + 7K*(y5 + 67 K*<y9 Z 
1 + 7K»<^ 2 + 87 2 K»</5 2 



- (l + 7K*^o) ^r, 



V 2 $ c 



-3Aa d 



1 + 7 (y*z<p + 
1 



B 2 H 2 \ 2 1 + 



1 + 71 \[kI<£ + 



+a 



•"0 -"0 



1 + 7 ( JK^if + "^ IL U 



1 + JK t ip 2 



7V K *<P 



1 + JK^ip 2 + 87 2 K*<p 2 



67 



77 2 + 77 2 



(BH ) 



-u + 3aAa 



1 



B 2 H?. 



(B3) 



(B4) 



where A is defined in Sect. 2.3 and is a constant of 0(1). written as 
We have also used the code unit for other quantities. 
This equation contains u, which must be solved from the 
scalar field equation of motion. 

The scalar field equation of motion can be similarly 
written. In order that the equation can be integrated 
into MLAPM, we need to discretise it for the application 
of Newton-Gauss-Seidel relaxation method. This means 
writing down a discrete version of this equation on a 
uniform grid with grid spacing h. Suppose we want to 
achieve second-order precision, as is in the default Pois- 
son solver of MLAPM, then V 2 it in one dimension can be in which 



V 2 u^V h2 Ui 



Uj+i + itj-i — 2uj 
h 2 



(B5) 



where a subscript j means that the quantity is evaluated 
on the j-ih point. The generalisation to three dimensions 
is straightforward. 

The discrete version of the equation of motion for u is 
then 



L h Kj, fc ) = 0, 



(B6) 



t h 1 \ ' 1 

L h (u i , j , k ) = - 



1 + •fK^lfi 2 + 67 2 K*<p 2 1 



+ 7K*<y5 2 h 2 

+ 3j^/K^lfi (l + 7K»(^g) fl 

+3a\a 3 



[ui+ij,k + u t-i,j,k + u>i,j+i,k + v>i,j-i,k + Ujj,fe+i + Ui^k-i - 6ui,j,k] 
Pc 1 



1 + 7 ( y/lU(p + -^-Ui,j,k 



1 + JK*lf 2 



1 



/ — I B2R l 



1 + Q 



, 1+a 



w + n 2 b 2 h 2 



67 



H 2 c 



2 iijik 



, . , , . B^m 

1 + 7 ( y/K^if + -T^-UiJ^ 



1 B2H ?> 



(1 + 3 7 ) 



-J7T- + 37V«*^— 772 



1 +7K*<^ (y/iuip) 



1 



/ B 2 H 2 V 

1 + 7 ( \f^ i P + —^-^uj.k \ 



1 + JK*(fi 2 



(B7) 



Then, the Newton-Gauss-Seidel iteration says that we 
can obtain a new (and often more accurate) solution of 



using our knowledge about the old (and less 



accurate) solution u° ld 



i,j,k 



old 

1 i,j,k 



via 



.old 



dL h 



.old 



/ du i,j,k 



(B8) 
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The old solution will be replaced by the new solution to 
Uij.k once the new solution is ready, using the red-black 
Gauss-Seidcl sweeping scheme. Note that 



d 



Ui,j, k 



1 + ~fK^ 2 + 6 7 V^ 2 6 U' + H 2 B 2 H 2 3a(l + a)\a 2 (BH /c) 

' 07 



1 + JK*<fi 2 



H 2 



( 1 , B2H l 



2+a 



B 2 H 2 

-67 2 ae2 V^*^ (l + 7 K *^o) flmPc 



l + 7 ( \J~k~*v + ^r-Uij^ 



B 2 H 2 
+12 7 2 ^-Xy^pa 3 



1 B 2 H* 
K*<P + -^r-U itji k 



ac* 



B 2 H 2 



+6aj 5— Xy/luipa" 



1 + 7 ( ^fnlp + 
1 



B 2 H% 



c 2 u i,j,k 



-47 



ac* 



2 B 2 H'q 



1 + 7 ( + tL ^r-u i j > k ) ( V^*<^ + -^r-Ui,j,k 



l+a 



(1 + 37)^35- + 3 7\/k^^-775 
^0 -"0 



B Hi 



K *f H ac 2 u i,j,k 



1 + 7 ( ^/kT^ + 



ac 2 



(B9) 



In principle, if we start from a high redshift, then the 
initial guess of Uij t k for the relaxation can be so chosen 
that the initial value of p in all space is equal to the 
background value p, because at this time we expect this 
to be approximately true any way. At subsequent time- 
steps we could use the solution for Uij k from the previous 
time-step as our initial guess. If the timestep is small 
enough then we would expect u to change only slightly 
between consecutive timesteps so that such a guess will 
be good enough for the iterations to converge quickly. 

ALGORITHM TO SOLVE THE BACKGROUND 
EVOLUTION 

Here we give our formulae and algorithm for the back- 
ground field equations which can also be applied to lin- 
ear Boltzmann codes such as CAMB. Throughout this Ap- 
pendix we use the conformal time r instead of the physi- 
cal time t, and ' = d/dr, % = a' /a. All quantities appear- 
ing here are background ones unless stated otherwise. 

For convenience, we will work with dimcnsionless quan- 
tities and define ip = y/i^Zp and N = In a so that 



(CI) 
(C2) 



* =n dN> 



4>"=u 



d 2 t/j 



dN 2 



dN' 



With these definitions it is straightforward to show that 
the scalar field equation of motion can be expressed as 



U_ 

Ho 



d 2 ip 



dN 2 
2 dV{i>) 



2% 
%o 



H 2 



dip 



1 a/2 + nj2 
Tin tin 



chp_ 
dN 

dip 



■0, (C3) 



where T-Ln is the current value of %. 

Obviously we need to know how to compute the quan- 
tities T-L/T-Ln and H' /Ho as well. For H/Ho, we start with 
the Fricdmann equation 



m 1 



1 



1 + / 
1 



[p m + Pr + V(lp)] a 
2 



+ 



1 + / 



dll> 

dN 



,df dtp 
'd^pdN 



•H 2 , (C4) 



where p m and p r are the energy densities for matter and 
radiation respectively. We define the fractional energy 
densities for matter and radiation respectively as 



_ KeffPmO 

3H 2 

K © OPrO 



1 



K*PmO 



1 + /0 



(C5) 



3H 2 



2(l + /o)+4(f)' 



K*PrO 



l + /0 



2(l + /o) + 3(4) 



3"Hq 



(C6) 



where a subscript means the present-day value. No- 
tice the difference between these definitions, which comes 
from the different treatments for radiation and matter in 
numerical codes such as CAMB. For radiation, e.g., pho- 
ton, we know the present temperature of the CMB and 
thus its exact energy density p r o, as well as the locally 
measured value of gravitational constant (which in 
scalar-tensor theories is in general different from k*) and 
current Hubble expansion rate Ho, and so the definition 
Eq. (C6) comes out naturally, where we have used the 
relation between k* and K0 1 . For matter, the fractional 
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energy density is to be interpreted from the cosmologi- 
cal observablcs such as CMB and large scale structure, 
which arc obviously different in ACDM and scalar-tensor 
theories; consequently there is some freedom in defining 
it and we make it as in Eq. (C5). 
Then, remembering that 



we have 



p m (xa 
p r oc a~ 



(C7) 
(C8) 



in which we have defined 



H 2 

Hn 



B: 



VL 

H 



ldf_H 2 

2 dip 



i i f i 1 df dip J at2 ' 
1 + J + Id^dN " JV 



(C14) 
(C15) 
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■£^n r a- 2 + (i + f ) n^r 1 + 

^0 



3%n 



iifl df dip 
* dp dN 



1 ( dip 
6 \ dN 



(C9) 



in which (where both k* and k^q are constants, and 
K0o is the present value of kq) 



2(l + /o) + 3(^ y 
(l + /o) ^ 



2(l + /o)+4($ 



(CIO) 



For H' /Ho, we use the Raychaudhrui equation 
1 



H' 



k* (p + 3p) a 2 

k* [p m + 2 p. 

d 2 f 



6 

1 1 

6TT7 

6TT7 ( 2 + 3 di 



2V(ip)] a 2 

d ^\ 2 nj 2 



1 1 



df ( 



2 1 + /dV V div 



(Cll) 



As in the above, dividing this by "H 2 and rearranging, we 
obtain 



H/_ 
H\ 



(i + fo) n^- 1 + 2^n r a~ 



1 _i_ f i \df_d±_ 

x ' J ' 2 dip dN 



1 df 

2 cijV 2 



2 ti^ 5 
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i , f , 1 d/ dip 
'J ~ 2 dip dN 



Ho 



(C12) 



Substituting Eqs. (C9, C12) into Eq. (C3), we finally 
arrive at 



1 



J ' 2 \dip 



1 4. f 4- idf_d± 
1 ' J ^ 2 dip dN 



A 



dN 2 



(2A + B] 



Kif dV 
~H%d^ ( 



dN 
df 



3(A + B)-f-=0, (C13) 



'dip 



where A and B do not contain d 2 ip/dN 2 , to lighten the 
notation. 

When solving for tp (or tp), we just use Eq. (C3) aided 
by Eqs. (C9, C12). It may appear then that, given any 
initial values for Vini and (dip/dN) ini , the evolution of tp 
is obtainable. However, Eq. (C9) is not necessarily sat- 
isfied for ip evolved in such way. Instead, it constrains 
the initial condition ip must start with, and the way it 
must subsequently evolve. This in turn is determined by 
the parameters A, a, 7, since a, 7 specify a model and are 
fixed once the model is chosen; the only concern is A. 

As for the initial conditions Vim and (dip/dN) ini , 
we have found that the subsequent evolution of V is 
rather insensitive to them. Thus, we choose Vim = 
(dip/dN) ini = at some very early time (say N mi cor- 
responds to o,j n i = e^'"' = 10~ 8 ) in all the models. Such 
a choice is clearly not only practical but also reasonable, 
given the fact that we expect that the scalar field starts 
high up the potential and rolls down subsequently. 

As for A, we use a trial- and-error method to find its 
value which ensures that (again subscript indicates the 
current time) 



-^-n r + (i + fo)n rn + 



1 + /o + (J0 o 



1 ( d±\ 

6 {dN) 



= 1 (C16) 



which comes from setting a = 1 in Eq. (C9). 

We determine the correct value of A for any given 
a, 7 in this way using a trial-and-error routine, and then 
compute the values of ip and dip/dN for predefined val- 
ues of N stored in an array. Their values at any time 
are then obtained using interpolation, and with these it 
is straightforward to compute other relevant quantities, 
such as H, H' and ip, which are used in the linear pertur- 
bation computations. 
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